home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qag.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-08-22  |  6.7 KB  |  266 lines
  1. /* integration/qag.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_errno.h>
  24. #include <gsl/gsl_integration.h>
  25.  
  26. #include "initialise.c"
  27. #include "set_initial.c"
  28. #include "qpsrt.c"
  29. #include "util.c"
  30.  
  31. static int
  32. qag (const gsl_function *f,
  33.      const double a, const double b,
  34.      const double epsabs, const double epsrel,
  35.      const size_t limit,
  36.      gsl_integration_workspace * workspace,
  37.      double * result, double * abserr,
  38.      gsl_integration_rule * q) ;
  39.  
  40. int
  41. gsl_integration_qag (const gsl_function *f,
  42.              double a, double b,
  43.              double epsabs, double epsrel, size_t limit,
  44.              int key,
  45.              gsl_integration_workspace * workspace,
  46.              double * result, double * abserr)
  47. {
  48.   int status ;
  49.   gsl_integration_rule * integration_rule = gsl_integration_qk15 ;
  50.  
  51.   if (key < GSL_INTEG_GAUSS15)
  52.     {
  53.       key = GSL_INTEG_GAUSS15 ;
  54.     } 
  55.   else if (key > GSL_INTEG_GAUSS61) 
  56.     {
  57.       key = GSL_INTEG_GAUSS61 ;
  58.     }
  59.  
  60.   switch (key) 
  61.     {
  62.     case GSL_INTEG_GAUSS15:
  63.       integration_rule = gsl_integration_qk15 ;
  64.       break ;
  65.     case GSL_INTEG_GAUSS21:
  66.       integration_rule = gsl_integration_qk21 ;
  67.       break ;
  68.     case GSL_INTEG_GAUSS31:
  69.       integration_rule = gsl_integration_qk31 ; 
  70.       break ;
  71.     case GSL_INTEG_GAUSS41:
  72.       integration_rule = gsl_integration_qk41 ;
  73.       break ;      
  74.     case GSL_INTEG_GAUSS51:
  75.       integration_rule = gsl_integration_qk51 ;
  76.       break ;      
  77.     case GSL_INTEG_GAUSS61:
  78.       integration_rule = gsl_integration_qk61 ;
  79.       break ;      
  80.     default:
  81.       GSL_ERROR("value of key does specify a known integration rule", 
  82.         GSL_EINVAL) ;
  83.     }
  84.  
  85.   status = qag (f, a, b, epsabs, epsrel, limit,
  86.                 workspace, 
  87.                 result, abserr, 
  88.                 integration_rule) ;
  89.   
  90.   return status ;
  91. }
  92.  
  93. static int
  94. qag (const gsl_function * f,
  95.      const double a, const double b,
  96.      const double epsabs, const double epsrel,
  97.      const size_t limit,
  98.      gsl_integration_workspace * workspace,
  99.      double *result, double *abserr,
  100.      gsl_integration_rule * q)
  101. {
  102.   double area, errsum;
  103.   double result0, abserr0, resabs0, resasc0;
  104.   double tolerance;
  105.   size_t iteration = 0;
  106.   int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
  107.  
  108.   double round_off;    
  109.  
  110.   /* Initialize results */
  111.  
  112.   initialise (workspace, a, b);
  113.  
  114.   *result = 0;
  115.   *abserr = 0;
  116.  
  117.   if (limit > workspace->limit)
  118.     {
  119.       GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
  120.     }
  121.  
  122.   if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
  123.     {
  124.       GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
  125.          GSL_EBADTOL);
  126.     }
  127.  
  128.   /* perform the first integration */
  129.  
  130.   q (f, a, b, &result0, &abserr0, &resabs0, &resasc0);
  131.  
  132.   set_initial_result (workspace, result0, abserr0);
  133.  
  134.   /* Test on accuracy */
  135.  
  136.   tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
  137.  
  138.   /* need IEEE rounding here to match original quadpack behavior */
  139.  
  140.   round_off = GSL_COERCE_DBL (50 * GSL_DBL_EPSILON * resabs0);
  141.  
  142.   if (abserr0 <= round_off && abserr0 > tolerance)
  143.     {
  144.       *result = result0;
  145.       *abserr = abserr0;
  146.  
  147.       GSL_ERROR ("cannot reach tolerance because of roundoff error "
  148.          "on first attempt", GSL_EROUND);
  149.     }
  150.   else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
  151.     {
  152.       *result = result0;
  153.       *abserr = abserr0;
  154.  
  155.       return GSL_SUCCESS;
  156.     }
  157.   else if (limit == 1)
  158.     {
  159.       *result = result0;
  160.       *abserr = abserr0;
  161.  
  162.       GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
  163.     }
  164.  
  165.   area = result0;
  166.   errsum = abserr0;
  167.  
  168.   iteration = 1;
  169.  
  170.   do
  171.     {
  172.       double a1, b1, a2, b2;
  173.       double a_i, b_i, r_i, e_i;
  174.       double area1 = 0, area2 = 0, area12 = 0;
  175.       double error1 = 0, error2 = 0, error12 = 0;
  176.       double resasc1, resasc2;
  177.       double resabs1, resabs2;
  178.  
  179.       /* Bisect the subinterval with the largest error estimate */
  180.  
  181.       retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
  182.  
  183.       a1 = a_i; 
  184.       b1 = 0.5 * (a_i + b_i);
  185.       a2 = b1;
  186.       b2 = b_i;
  187.  
  188.       q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
  189.       q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);
  190.  
  191.       area12 = area1 + area2;
  192.       error12 = error1 + error2;
  193.  
  194.       errsum += (error12 - e_i);
  195.       area += area12 - r_i;
  196.  
  197.       if (resasc1 != error1 && resasc2 != error2)
  198.     {
  199.       double delta = r_i - area12;
  200.  
  201.       if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
  202.         {
  203.           roundoff_type1++;
  204.         }
  205.       if (iteration >= 10 && error12 > e_i)
  206.         {
  207.           roundoff_type2++;
  208.         }
  209.     }
  210.  
  211.       tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
  212.  
  213.       if (errsum > tolerance)
  214.     {
  215.       if (roundoff_type1 >= 6 || roundoff_type2 >= 20)
  216.         {
  217.           error_type = 2;    /* round off error */
  218.         }
  219.  
  220.       /* set error flag in the case of bad integrand behaviour at
  221.          a point of the integration range */
  222.  
  223.       if (subinterval_too_small (a1, a2, b2))
  224.         {
  225.           error_type = 3;
  226.         }
  227.     }
  228.  
  229.       update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
  230.  
  231.       retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
  232.  
  233.       iteration++;
  234.  
  235.     }
  236.   while (iteration < limit && !error_type && errsum > tolerance);
  237.  
  238.   *result = sum_results (workspace);
  239.   *abserr = errsum;
  240.  
  241.   if (errsum <= tolerance)
  242.     {
  243.       return GSL_SUCCESS;
  244.     }
  245.   else if (error_type == 2)
  246.     {
  247.       GSL_ERROR ("roundoff error prevents tolerance from being achieved",
  248.          GSL_EROUND);
  249.     }
  250.   else if (error_type == 3)
  251.     {
  252.       GSL_ERROR ("bad integrand behavior found in the integration interval",
  253.          GSL_ESING);
  254.     }
  255.   else if (iteration == limit)
  256.     {
  257.       GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);
  258.     }
  259.   else
  260.     {
  261.       GSL_ERROR ("could not integrate function", GSL_EFAILED);
  262.     }
  263. }
  264.  
  265.  
  266.